Gravitational waves from oscillating accretion tori: 
Comparison between different approaches 
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Quasi-periodic oscillations of high density thick accretion disks orbiting a Schwarzschild black 
hole have been recently addressed as interesting sources of gravitational waves. The aim of this 
paper is to compare the gravitational waveforms emitted from these sources when computed using 
(variations of) the standard quadrupole formula and gauge-invariant metric perturbation theory. To 
this goal we evolve representative disk models using an existing general relativistic hydrodynamics 
code which has been previously employed in investigations of such astrophysical systems. Two are 
the main results of this work: First, for stable and marginally stable disks, no excitation of the 
black hole quasi-normal modes is found. Secondly, we provide a simple, relativistic modification of 
the Newtonian quadrupole formula which, in certain regimes, yields excellent agreement with the 
perturbative approach. This holds true as long as back-scattering of GWs is negligible. Otherwise, 
any functional form of the quadrupole formula yields systematic errors ~ 10%. 

PACS numbers: 04.30.Db, 95.30.Lz, 98.62.Mw 
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Quasi-periodic oscillations of axisymmetric thick ac- 
cretion disks or tori orbiting around a black hole have 
been recently addressed as promising sources of gravita- 
tional waves (GWs) Q. In Ref. Q the GW signals were 
computed by means of the quadrupole stress formula 
within the Newtonian approximation, as it is commonly 
done in the literature 2, 3, 0,1^1^. The purpose of the 
present article is to compare the gravitational waveforms 
computed using the quadrupole formalism with those 
obtained from the theory of (gauge-invariant) gravita- 
tional metric perturbations of Schwarzschild spacetime, 
which amounts to solve (in the time domain) the inhomo- 
geneous even- (Zerilli-Moncrief) and odd-parity (Regge- 
Wheeler) master equations 0, H, • Since our interest is 
on the comparison between the two approaches, we shall 
concentrate on even-parity perturbations only. We note 
that similar investigations have been reported by [g in 
the context of oscillations of neutron stars, and by [W^ 
for orbiting particles around a Schwarzschild black hole. 



three odd-parity multipoles 
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with the even multipoles transforming as (—1)' and the 
odd as (—1)'^^ under a parity transformation {9,ip) — > 
{■n — 6, 77 + If). The presence of a fluid matter disk around 
the black hole is accounted through a "source" term 
in the linearized Einstein's equations, represented by a 
stress energy tensor t„,p which can be expanded in mul- 
tipoles as well (see [ll|, [l^ for details) . By taking suit- 
able combinations of the metric multipoles, it is possible 
to define an even-parity Z^™ (Zerilli-Moncrief) and an 
odd-parity 1?''^™ (Regge- Wheeler) gauge-invariant mas- 
ter function. The linearized Einstein's equations for the 
even-parity part reduce to the Zerilli-Moncrief equation 
in Z'^" 
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In order to perform such a comparison we use the same 
physical setup of [H , namely a relativistic thick accretion 
disk undergoing oscillations induced by suitable initial 
perturbations of an equilibrium configuration. A num- 
ber of restrictive hypothesis are assumed in our compu- 
tational framework which, however, are valid as long as 
the mass of the disk fi is much smaller than the mass 
of the black hole AT. In such a case, the self-gravity 
of the disk and the radiation reaction effects can be ne- 
glected. We assume thus that the spacetime metric g^i, 
is described by the Schwarzschild metric g^i, plus a non- 
spherical perturbation h^„. As a result of the spherical 
symmetry, /i^^ can be expanded in seven even-parity and 



where r* = r + 2M \og[r / (2M) — 1] is the Regge- Wheeler 
tortoise coordinate and 
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A(A - 2)^r-^ + 6(A - 2)^MH + 36(A - 2)M^r + 72AP 
r3 [(A - 2)r + 6M]^ 

is the Zerilli potential, with A = + Explicit expres- 
sion for the source 5*1™, in the gauge- invariant formalism 
of Gerlach-Sengupta JJj, can be found in Refs. 11, 12j. 
In the wave-zone (i.e. for large r) it can be shown that 
the even-parity GW amplitude reads (see e.g. ^2] and 
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TABLE I: Marginally stable (A, C) and stable (B, D) constant angular momentum thick disks orbiting around a Schwarzschild 
black hole of mass M = 2.5Mq. From left to right the columns report the name of the model, the number of radial and polar 
grid zones used in the simulations, the disk-to-hole mass ratio, the polytropic constant k of the isoentropic equation of state 
(EOS) p = Kp'' with 7 — 4/3, the value of the specific angular momentum the position of the cusp rcusp and of the center 
?"centor of the disk, the density at the center pc, the location of the inner (rin) and outer (rout) disk boundaries, the linear size 
in the equatorial plane L, the orbital period of the center torb, the value of the potential barrier AW, and the magnitude of 
the perturbation parameter r/. The inner and outer radii of the computational domain are rmin — 2.1 (for all models) and 
rmax = 18 (A), 20 (B and C), and 40 (D). All radii are given in units of M. 



Model 


Nr 


Ne 


^^/M 


K (cgs) 


I 


^cusp 


^center 


Pc (cgs) 


Tin 


rout 


L (km) 


torb (ms) 


AW 


V 


A 


320 


150 


0.10 


3.01 X lO" 


3.75 


4.835 


7.720 


3.58 X 10" 


4.83 


11.57 


24.87 


1.66 





0.1 


B 


300 


150 


0.02 


3.58 X 10^3 


3.75 


4.835 


7.720 


1.39 X 10" 


5.81 


10.32 


16.65 


1.66 


-0.002 


0.01 


C 


300 


150 


0.10 


9.60 X 10" 


3.80 


4.576 


8.352 


1.15 X 10" 


4.57 


15.89 


41.76 


1.86 





0.01 


D 


300 


150 


0.01 


6.19 X 10" 


3.95 


4.107 


9.971 


2.83 X 10" 


5.49 


29.08 


87.09 


2.43 


-0.015 


0.02 



references therein) 

1=2 m=-l V ^ 
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where _2l^^'" (^i v) ^^'^ the s = — 2 spin- weighted spher- 
ical harmonics. In axisymmetry {m = 0) 

1=2 

where = Yjg% - cot^F^o, and /ix = 0. The total 
energy radiated in GWs is defined by 

where we have introduced the Fourier transform 
Z{lj, r) = e-"^*Z(t, r)dt, with uj = 2nf and Z = 
In the following we focus our discussion on the / = 2 
multipole where = A/45/(47r) sin^ 9. 

The equilibrium accretion disks that we use as initial 
models for our simulations, and whose hydrodynamics is 
accounted for in the source terms of Eq. , are built by 
solving the relativistic Bernoulli equation as described 
in We only recall here that the tori are isoentropic, 
non-Keplerian (thick) disks with a cusp on the equato- 
rial plane, through which matter can accrete onto the 
black hole provided there is some departure from hy- 
drostatic equilibrium. For simplicity, we only consider 
initial models with a uniform radial distribution of the 
specific angular momentum, a fact which guarantees that 
the radial location of the cusp is always smaller than the 
marginally stable orbit. Despite constant specific angu- 
lar momentum disks do not represent physical tori, be- 
ing dynamically unstable against nonaxisymmetric insta- 
bilities, such a simple model is still appropriate for our 



purpose of testing the quadrupolc formula (see [lj| for 
details on improved models of equilibrium tori) . The ge- 
ometrical structure of the equipotential surfaces of the 
tori resembles that of the Roche lobes of a binary sys- 
tem. In particular, by varying the value of the potential 
barrier AW = Win — W^cusp at the inner edge of the disk, 
it is possible to build either a stable configuration corre- 
sponding to a torus inside its Roche lobe (AW^ < 0), or 
a torus overflowing its Roche lobe (AW^ > 0) and there- 
fore accreting onto the black hole. The case AW — 
corresponds to the marginally stable case. 

In the present work we deal with models which are 
stable or marginally stable with respect to the runaway 
instability (see e.g. and references therein). Table HI 
lists the main properties of the models considered. The 
mass of the black hole is M = 2.5Mq, and disks with dif- 
ferent disk-to-hole mass ratio fi/M are considered. In all 
cases, but particularly for models B and D, the disk rep- 
resents a small perturbation to the gravitational field of 
the black hole. Models A and B share the same location 
of the center and the cusp (i.e. have the same equipoten- 
tial surfaces), but differ in the total mass /x and in the 
value of the potential gap AW. In order to induce peri- 
odic oscillations of the disks we follow the same procedure 
of P] and introduce a small radial velocity perturbation, 
parameterizing its strength in the form Vr — rjv^P^, where 
w^P^ is the (stationary) radial velocity of the spherically 
symmetric low density atmosphere surrounding the torus. 
The values of the perturbation parameter chosen are typ- 
ically in the interval r/ e [0.01,0.1]. This produces low 
mass accretion rates (see below) and validates the as- 
sumption of a "fixed" background spacetime. The initial 
data are evolved in time using an axisymmetric general 
relativistic hydrodynamics code based on the conserva- 
tive formulation of the equations, whose accuracy has 
been assessed in earlier investigations Umill. The 
computational (r, 0)-grid used in the simulations has non- 
uniform radial spacing and equidistant polar spacing. Ta- 
ble ^ reports the number of grid zones employed for each 
model. 

It needs to be mentioned that setting up the initial data 
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FIG. 1: Left panel: Time evolution of the I = 2,m = GWs for models A and B, computed using the Zerilli-Moncrief equation. 
The first 60 periods of a t = lOOtorb simulation are shown. The oscillatory pattern in the waveforms is the response of the 
quasi-periodic oscillations of the disks. Right panel: Energy spectra for model A (solid line) and B (dashed line). Only fluid 
p-modes are visible in the spectra (enlarged in the inset, where a Hamming filter was used) which show no evidence of the 
excitation of the fundamental QNM of the black hole (its frequency is indicated by a vertical dotted line at 2Muj — 0.7473). 



for even-parity perturbations requires some care [ll|, as 
the initial profile of Z and Z,f should be obtained by 
solving the (linearized) Hamiltonian and momentum con- 
straints with matter sources, respectively. When this is 
not done, an initial unphysical burst of gravitational radi- 
ation is produced, which may cause the excitation of the 
black hole quasi-normal modes (QNM) and hide physical 
effects if these were to happen on a comparable timescale. 
However, this is not the case in our simulations as the 
timescale of the oscillations of the disk is much larger 
than the decay time of the black hole QNM. Hence, once 
the initial data have been evolved for a sufficient time 
(roughly a couple of orbital periods) the gravitational 
waveforms remain entirely unaffected by the initial GW 
content. In practice, we choose for the initial profile of Z 
a. I — 2 time-symmetric stationary multipole by solving 
Eq. lO as an ODE with Z^tt = 0. 

The left panel of Fig. ^ shows the I — 2 gravitational 
waveforms extracted for models A and B through the 
Zerilli-Moncrief equation (at robs = 400M) as a function 
of the observer retarded time u = t — r°^^, in units of the 
orbital period of the center of the disk. Due to the fi- 
nite location of the observer a vertical offset is present in 
the waveforms, induced by the matter sources in Eq. |(2Jl. 
This offset has been corrected for in the figure. The to- 
tal mass accreted onto the black hole at the end of the 
simulations (t = lOOtoib) is ~ 0.03M for model A and 
^ 2.20 X lO^^A/ for model B, which justifies the fixed 
spacetime approximation we use. The waveforms shown 
in the left panel of Fig. ^ indicate that the GW signals 
follow the torus oscillations. The decrease of the GW am- 
plitude is determined by the amount of accreted mass. 
For model B the GW amplitude remains almost con- 
stant throughout the evolution, while for model A, where 



more mass transfer is possible at the cusp, it decreases 
by ^ 60% after 100 orbital periods. As becomes evident 
on the right panel of the figure, which shows the corre- 
sponding energy spectra for both models, no black hole 
QNM is excited in the process. The spectra have been 
obtained after cutting out the first four periods of oscil- 
lation in order to avoid any contamination from the early 
time unphysical burst. At low frequencies (see the inset 
in the right panel of Fig.QJ, distinctive peaks are visible, 
which correspond to the fluid p-modes of the disk, namely 
a fundamental mode / and a series of overtones P, [Tg . 
The fundamental I — 2 QNM of the black hole should lie 
at 2Mllj = 0.7473 (indicated by a vertical dotted line in 
the plot), but the computed spectra show no hints of this 
mode. This is consistent with the results reported in 
for the case of radial infall of quadrupolar shells, where it 
was shown that the fundamental QNM of the black hole 
is only significantly excited if the bulk of the maximum 
density of the (compact) shell reaches the peak of the 
Zerilli potential and if its typical linear dimensions are 
comparable with the width of the potential itself. 

The absence of black hole QNM excitation in our simu- 
lations suggests that the gravitational waveforms should 
also be accurately computed by means of the Landau- 
Lifshitz quadrupole formula. We use our setup to assess 
the validity of the quadrupole formalism through a di- 
rect comparison with the gauge-invariant GWs extrac- 
tion shown in Fig. ^ Since the quadrupole formula has 
been derived in a Newtonian framework, a number of 
approximations must be first considered in the transi- 
tion to a general relativistic context. Let us consider a 
Minkowski background spacetime r]f^i, (/t — 0, 1, 2, 3) and 
an axisymmetric GW source whose physical dimensions 
are much smaller than the reduced gravitational wave- 
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ordinates we obtain 




FIG. 2: Comparison between the I = 2 gravitational wave- 
forms for model B extracted using the quadrupole formula 
SQFi (dotted line) and the gauge-invariant approach (solid 
line). The bottom panel shows the relative difference be- 
tween the two normalized to the mean amplitude of the gauge- 
invariant method. 



length (i.e. slow- motion approximation). In such case the 
only non- vanishing polarization amplitude of the wave is 
which is given by 6, 16] 



_ Ssin^e d^j:.. 



where 



(7) 



(8) 



becomes the reduced matter quadrupole moment tensor 
in the Newtonian limit = p , p being the rest-mass 
density of the source. Eq. ((Tj) is known as the Standard 
Quadrupole Formula (SQF) and can be used to extract 
GWs in a slow-motion nearly Newtonian regime. An ex- 
tension of the simple Eqs. (TJ and © to general relativity 
is problematic because the Einstein's equations are non- 
linear in the background metric g^j/ 16j. However, since 
in the present physical context no black hole QNMs ex- 
citation is found, the nonlinear structure of the back- 
ground metric, represented by the potential can 
be neglected in practice and the use of the quadrupole 
formula can be justified a posteriori. Following a prag- 
matic approach, a reasonable modification to account for 
general relativistic effects can be obtained by projecting 
i^" onto the local inertial frame of the tetrad associated 
with the stationary observers of the Schwarzschild space- 
time, so that = a^t°°, where t' 
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and a ~ \J~-gm is the lapse. Here p is the pressure, 
h = 1 + e+ p/pis the specific enthalpy, and e is the 
specific internal energy of the fluid. Using spherical co- 



/ = 27ry dz [hDW -v){^z^ -^r^dr 



(9) 



where / = 8^33/2, z — cos9, D = pW is the relativistic 

(1 - gijv'v')-^^^ is 



au 



rest-mass density, and W 
the Lorentz factor. If we define Af(f = y^MnJTE d'^I/dt^, 
the energy emitted in the quadrupole approximation is 



E 



327r2 



u'\A^S{Lo,r)\'duj 



(10) 



In the following we use SQFi to denote the standard 
quadrupole formula with the generalized quadrupole mo- 
ment proposed above. Figure [3 displays the comparison, 
for model B, of the gravitational waveforms obtained 
with the Zerilli-Moncrief equation and with the SQFi 
definition for J- 33. The second time derivative has been 
computed directly from the numerical data using cen- 
tered finite-differences. We note that for our particular 
physical setup and unlike the case of core collapse simula- 
tions 0, J the numerical noise associated with this pro- 
cedure is negligible. Fig. |21 shows an excellent agreement 
between the two methods: the two lines perfectly overlap 
at early times and only small differences are reached af- 
ter many orbital periods. Fig.^lshows the corresponding 
comparison of the energy spectra for t = lOOiorb and (in 
the inset) a close up of the waveforms for different alter- 
native definitions of The Newtonian SQF (i°° = p) 
slightly underestimates the wave amplitude while SQF2 
(t°" = pW) is closer to the gauge-invariant signal; SQF3, 
corresponding to = pvP and used in Refs. 0,|E0|j 
gives the worst relative agreement in the present physical 
context (dash-dotted line in the inset of Fig. |3Jl. 

To quantify the differences among the various alterna- 
tive expressions we use for the quadrupole moment we 
report in Table ITU the relative errors Ai?^*^ found in the 
total energy emitted after the frequency integration of 
the power spectra for all four models of our sample. We 
recall that from model A to D the radial location of the 
center of the disk increases gradually. It is found that the 
SQFi (fourth column in Table ITl|l gives the best agree- 
ment with the gauge-invariant GW extraction, with an 
error that is always less than 2%. On the other hand, the 
SQF3 yields the worst result for either model, systemati- 
cally overestimating the amplitude of the waveforms, the 
reason being that, when compared to the other expres- 
sions, the functional form of SQF3 contains a division by 
an extra lapse function (a < 1). 

It needs to be mentioned, however, that the decrease in 
Ai5^° with increasing reenter, observed in SQF and SQF3 
in Table seems accidental. This has been checked 
by analyzing additional models with the same linear di- 
mension but with the center located at reenter = 10.17 
and reenter = 10.47. The latter is the maximum al- 
lowed value for a closed disk with I = 4 No clear 
continuation in the trend in Ai?^" reported in Table ^ 
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FIG. 3: Comparison of the energy spectra for model B and 
a portion of the waveform when different definitions of the 
generaUzed quadrupole moment of the torus are used. 




FIG. 4: Model B: comparison between the Zerilli-Moncrief 
function, the SQFi and the SF for different choices of the 
potential The inset shows more periods of oscillations. 



TABLE II: Relative differences AS^" = \E'^° ~ E\/ E'^° a.raong, 
the emitted energies computed using different versions of the 
quadrupole formula {E) and the gauge-invariant approach 
(i?'^") for disk models with different position of the center. 



Model 


"^center 










A 


7.720 


2.9% 


0.7% 


0.8% 


11.0% 


B 


7.720 


5.9% 


0.3% 


2.6% 


10.6% 


C 


8.352 


4.9% 


1.8% 


3.2% 


6.7% 


D 


9.971 


2.7% 


1.7% 


2.2% 


3.8% 



has been found for these additional models. The largest 
errors computed appear in SQF3 for the model with 
^center = 10.17 (AS^o - 5.4%) and in SQF for the model 
with reenter = 10.47 {AE^° - 5%). 

Further evidence on the sensitivity of these results on 
the model properties is gained by reducing the potential 
gap in model D to AW = —0.008, keeping the center at 
the same location. As a result the disk becomes much 
larger, extending from Hn = 4.93 to rout =42.5. In this 
case, the relative difference in the emitted energy be- 
tween the gauge-invariant waveform and the quadrupole 
formula can be as large as ~ 25% for SQF, SQFi and 
SQF2, while A£;20 - 21% for SQF3. Following Tanaka 
et al. 10], who computed GWs from particles orbiting 
around a Schwarzschild black hole using linear perturba- 
tion theory, we interpret this disagreement for enlarged 
tori as the result of GWs back-scattering with the tail of 
the background potential. 

As mentioned above, in simulations of core collapse 0, 
1^, the accurate numerical computation of the time 
derivatives appearing in the SQF is prone to introduce 
high-frequency numerical noise. The common procedure 
to circumvent this is to rewrite Eq. Q making use of 



the Newtonian limit of the general relativistic hydrody- 
namics equations (i.e. the Euler and continuity equa- 
tions) 0, [13 to eliminate the time derivatives. In ax- 
isymmetry the resulting expression reads 



A 



E2 



k I dz J p 



VrV''i3z^-l) 



r^dr , 



r9^$(3z^-l)+ 3z<9e$Vl-22 



(11) 



where k — IQir^/"^ / and <& is the Newtonian gravita- 
tional potential generated by the central source M . This 
equation is known as the stress formula (SF), and it has 
been applied in diverse astrophysical scenarios 0, 0, . 
Since Eq. is obtained from the Newtonian fluid equa- 
tions, its application in a general relativistic code poses 
some ambiguities. However, it is possible to introduce 
suitable forms for $ in Eq. (|ll|l to obtain results which 
are consistent with those derived in the perturbative 
framework. An elementary possibility is to use a first 
post-Newtonian (IPN) expression for $, as it was done 
in 1]. In practice <I> can be obtained from g^r = 1 — 2$ 
which is correct at IPN [3; i.e. $ = = -M/(r - 

2Af). Moreover, one can also introduce the first terms of 
the 2PN expansion of Qrr ^3 ' solving the quadratic equa- 
tion Qrr = 1 — 2$ -I- 2$2 and choosing the solution which 
recovers the Newtonian potential $^ = —M/r for large r; 
i.e, $ = $2PN ^ 1/2 [1 - [r + 2M) / {r - 2M)\ . Fig- 
ure^ shows the GW signals for model B computed using 
the SQFi (long-dashed line), the Zerilli-Moncrief func- 
tion Z (sohd Hue), and the SF with (dotted Hue) 
or $2PN (ghort-dashed line). The waveforms obtained 
from the SF present a systematic global vertical offset 
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which was corrected in the figure. Figure 0] shows that a 
good agreement is found among all waveforms except for 
those computed using In this case, the signal is 

out of phase by half a period, due to the dr^ term being 
dominant in Eq. Hill) and thus responsible for a change of 
sign, while the amplitude is about a factor of two larger 
than that of the other waveforms. On the other hand, 
it is remarkable the good agreement in amplitude found 
among the SF waveform using with both the Zerilli- 
Moncrief equation and the SQFi. We finally note that 
the use of a potential consistent with the Schwarzschild 
metric, obtained from goo — exp(— 2$), results in a signal 
whose amplitude is ~ 30% larger than the one obtained 
with SQFi. This stresses the intrinsic ambiguity which 
relies in the use of a modified Newtonian formula in a 
general relativistic context. Work is in progress to derive 
a modified stress formula which holds in this regime. 



In this paper we have studied the GW signals 
driven by oscillations of thick accretion disks around a 
Schwarzschild black hole by means of gauge-invariant 
metric perturbation theory and modifications of the SQF. 
The comparison between the two approaches used has 
been discussed in detail. For the models considered, the 
GW signal has been found to be mainly determined by 
the matter fiows and no evidence of black hole QNM ex- 
citation has been observed. 
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